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ABSTRACT 

The current picture of terrestrial planet formation relies heavily on our understanding of the dynam- 
ical evolution of planetesimals — asteroid-like bodies thought to be planetary building blocks. In 
this study we investigate the growth of eccentricities and inclinations of planetesimals in spatially 
homogeneous protoplanetary disks using methods of kinetic theory. Emphasis is put on clarifying the 
effect of gravitational scattering between planetesimals on the evolution of their random velocities. 
We explore disks with a realistic mass spectrum of planetesimals evolving in time, similar to that 
obtained in self-consistent simulations of planetesimal coagulation: distribution scales as a power law 
of mass for small planetesimals and is supplemented by an extended tail of bodies at large masses 
representing the ongoing runaway growth in the system. We calculate the behavior of planetesimal 
random velocities as a function of the planetesimal mass spectrum both analytically and numerically; 
results obtained by the two approaches agree quite well. Scaling of random velocity with mass can 
always be represented as a combination of power laws corresponding to different velocity regimes 
(shear- or dispersion-dominated) of planetesimal gravitational interactions. For different mass spectra 
we calculate analytically the exponents and time dependent normalizations of these power laws, as 
well as the positions of the transition regions between different regimes. It is shown that random 
energy equipartition between different planetesimals can only be achieved in disks with very steep 
mass distributions (differential surface number density of planetesimals falling off steeper than m~ 4 ), 
or in the runaway tails. In systems with shallow mass spectra (shallower than m -3 ) random velocities 
of small planetesimals turn out to be independent of their masses. We also discuss the damping ef- 
fects of inelastic collisions between planetesimals and of gas drag, and their importance in modifying 
planetesimal random velocities. 

Subject headings: planetary systems: formation — solar system: formation — Kuiper Belt 



1. INTRODUCTION. 

Gravity-assisted agglomeration of planetesimals in the 
protoplanetary disks around young stars is a backbone 
of a currently widely accepted paradigm of terrestrial 
planet formation, despite many uncertainties in the de- 
tails of specific processes involved. In many respects our 
understanding of the planetesimal disk evolution heavily 
relies on the numerical investigations, and these studies 
significantly grew in complexity in recent years (Wetherill 
& Stewart 1989, 1993; Kenyon & Luu 1998, 2002; Inaba 
et al. 2001). On one hand, the sheer scale of simulated 
systems has expanded considerably, which was made pos- 
sible by the increase in computing power; on the other 
hand, a wider range of concomitant physical phenomena 
(fragmentation, migration, etc.) can now be routinely 
followed simultaneously with the evolution of planetesi- 
mal mass and velocity distributions. 

Unfortunately, very often this complexity creates a 
problem when one tries to understand a relative role of 
each of the specific processes in shaping the planetesimal 
mass and velocity spectra. Disentangling the contribu- 
tions of different physical nature in the results of numer- 
ical simulations certainly is a tantalizing task which can 
often yield very confusing conclusions. Thus it is very im- 
portant to be able to isolate different physical processes 
operating in the planetesimal disks and to study their ef- 
fects on the disk evolution separately. One way to do this 
is to build simple but realistic models which are easy to 
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analyze and at the same time are able to grasp the main 
features of the physics involved. Such approach would 
grant us good analytical understanding of relevant pro- 
cesses, and we are going to follow this route in our study. 

Behavior of planetesimal velocities is one of the cru- 
cial ingredients of the terrestrial planet formation pic- 
ture because of the very important role played by the 
gravitational focusing (which depends on planetesimal 
velocities) in assembly of the big bodies in planetesimal 
disks. The purpose of this paper is to explore the evo- 
lution of inclinations and eccentricities (which represent 
planetesimal random velocities) caused by the mutual 
gravitational scattering of planetesimals. This is a clean 
and well-posed problem. The effects of gas drag and in- 
elastic collisions are initially neglected but later on we 
discuss how their inclusion affects our results. To keep 
things as clear as possible we have opted not to try cou- 
pling self-consistently planetesimal velocity evolution to 
the evolution of the mass spectrum in this study. Instead, 
we follow changes of random velocities in a disk with a 
distribution of planetesimal sizes which evolves in some 
prescribed manner. However simplistic this assumption 
seems to be, in adopting it we try to use mass spectra 
which are similar to those produced by self-consistent 
coagulation simulations, and to study a wide variety of 
such mass distribution models. 

Another argument in favor of this approach is the fact 
that the dynamical relaxation time in a disk of gravita- 
tionally interacting planetesimals is much shorter than 
the physical collision timescale (characterizing evolution 
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of the mass spectrum) if relative random velocities of 
planetesimals are smaller than the escape velocities from 
their surfaces (which is very often the case) . The dynam- 
ical evolution of planetesimal disk can then be approxi- 
mated as that of a system with a quasi-stationary plan- 
etesimal mass spectrum. Understanding this problem is 
a logical step necessary to provide a clearer perspective 
on how to build a fully self-consistent theory of planet 
formation. 

We describe the statistical approach to the problem 
of the evolution of planetesimal random velocities in [JU 
Our model of planetesimal mass spectrum is introduced 
In 21 we outline the results for the distribution of 
planetesimal random velocities vs. planetesimal masses 
as functions of the input mass distribution and time; we 
also compare these analytical predictions with numer- 
ical results. Table IC1I summarizes our major findings 
in a concise form. In S|2]we comment on how gas drag 
and inelastic collisions between planetesimals can affect 
planetesimal velocity spectra. We conclude in JS]by the 
discussion of our results and comparison with previous 
studies. 

2. VELOCITY EVOLUTION EQUATIONS. 

When the number of planetesimals in the protoplane- 
tary disk is very large and their masses are not sufficient 
to induce strong local nonuniformities in the disk (Ida 
& Makino 1993; Rafikov 2001, 2003b, 2003c), statistical 
approach and homogeneous "particle in a box" assump- 
tion are very helpful in the treatment of planetesimal 
evolution (Wetherill & Stewart 1989). We assume that 
planetesimals with mass to have velocity dispersions of 
eccentricity and inclination a e (m,t) and cr^(m, t) associ- 
ated with them. It is natural to characterize planetes- 
imal disk by its mass distribution, and we assume that 
N(m, t)dm is the surface number density of planetesimals 
with masses between to and to + dm. 

In studying the dynamics of planetesimal disks we 
make the following natural assumptions: 

• disk is locally homogeneous in radial direction and 
azimuthally symmetric; it is in Keplerian rotation 
around central star, 

• epicyclic excursions of planetesimals in horizontal 
and vertical directions are small compared with the 
distance to the central object, 

• planetesimal masses are much smaller than the 
mass of the central object M c . 

Two important consequences immediately follow from 
this set of assumptions. First, the gravitational scatter- 
ing of planetesimals during their encounter can be stud- 
ied in the framework of Hill approximation (Henon & 
Petit 1986; Ida 1990; Rafikov 2003a). In this approach, 
gravitational interaction between two bodies with masses 
to and to* introduces a natural lengthscale — Hill ra- 
dius 1 
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1 This definition follows original work of Henon & Petit (1986) 
and differs from Rjj = a[(m + m*)/3Af c ] 1 / 3 often used in the 
literature (e.g. Ida 1990; Stewart and Ida 2000). 
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where a is the distance from the central object. Nu- 
merical estimate made in Q assumes M c = Mq and 
<Xau = AU), and is intended to illustrate that typi- 
cally Rh <C a. Hill approximation yields two significant 
simplifications: 

• The outcome of the interactions between two bod- 
ies depends only on their relative velocities and dis- 
tances. 

• If all relative distances are scaled by Rjj and rel- 
ative velocities of interacting bodies are scaled by 
QRh, then the outcome of the gravitational scat- 
tering depends only on the initial values of scaled 
relative quantities. 

Second, numerical studies (Greenzweig & Lissauer 
1992; Ida & Makino 1992) have demonstrated that in ho- 
mogeneous planetesimal disks the distribution function 
tp(e,i) of absolute values of planetesimal eccentricities 
e and inclinations i is well represented by the Rayleigh 
distribution: 
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where er e and <Ji are the aforementioned dispersions of ec- 
centricity and inclination. It also follows from azimuthal 
symmetry that horizontal and vertical epicyclic phases r 
and u are distributed uniformly in the interval (0, 2it). 
These facts have important ramifications as we demon- 
strate below. 

Let's consider the gravitational interaction of two plan- 
etesimal populations: one with mass m and eccentricity 
and inclination dispersions a e and oi and the other with 
mass to* and dispersions <r* and a*. When the veloc- 
ity 2 distribution function of planetesimals has the form 
(J2J) it can be shown (Rafikov 2003a) that the evolution 
of cr e and <Tj due to the gravitational interaction with 
population of mass to* proceeds according to 
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where N* is the surface number density of planetesimals 
of mass to*, and A = —(r/2)d£l/dr is a measure of shear 
in the disk [in Keplerian disks A — —(3/4)0]. Scattering 
coefficients H± t 2 are defined as 

Hi = Hi (a er ,a ir ) 



de r di r ^ r (e r ,i r ) / dh\h\ (Ae 



2 Throughout the paper we often refer to eccentricities and in- 
clinations of planetesimals as their random velocities. 
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H 2 = H 2 (a- er ,<7 ir ) 

de r di r ip r (e r ,i r ) I dh\h\(e r • Ae sc ), (4) 



where a er ,ir = {o'ei + (T * 2 i) 1 ^ 2a l Rh are the dispersions 
of relative eccentricity and inclination normalized in Hill 
coordinates [see JjQl], h is a semimajor axes difference 
of interacting bodies scaled by Rh, and e r ,i r are the 
relative eccentricity and inclination vectors in Hill units 
(Goldreich & Tremaine 1980; Ida 1990). Function Ae sc 
represents a change of e r produced in the course of scat- 
tering, which is a function of not only absolute values of 
e r , i r , and h, but also of the relative epicyclic phases (for 
which reason we use vector eccentricities and inclinations 
here). The distribution function of relative eccentricities 
and inclinations of planetesimals t/j r {e r ,i r ) is analogous 
in its functional form to @, but with a e ^ replaced by 
relative velocity dispersions a er ^ r . Expressions analogous 
to EJ} can be written down also for inclination scattering 
coefficients K\ and K 2 . Note that the only assumption 
used in deriving is that of the specific form ijSJ) of the 
distribution function of planetesimal eccentricities and 
inclinations (Rafikov 2003a). 

Equations © describe dynamical evolution driven only 
by gravitational scattering of planetesimals, implicitly 
assuming that they are point masses. Physical size of 
planetesimal r p « 3 x 10~ 7 AU (m/(10 21 g)) 1 / 3 (for 
physical density 3 g cm -3 ) is very small compared to 
the corresponding Hill radius JJ}, which often justifies 
the neglect of physical collisions between planetesimals. 
However, at high relative velocities — higher than the 
escape speed from planetesimal surfaces — one can no 
longer disregard highly inelastic physical collisions which 
strongly damp planetesimal random motions. In this 
study we proceed by assuming first that velocities of 
planetesimals are below their escape velocities, and then 
discussing in JO how abandoning this assumption affects 
our conclusions. 

Equation Q and definitions (@J have been previously 
derived by other authors in a slightly different form 
(Wetherill & Stewart 1988; Ida 1990; Tanaka & Ida 1996; 
Stewart & Ida 2000): 

4 / 3 ™* 
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m + nV <7g + crj^ 
In this form the first term in brackets on the r.h.s. de- 
scribes the so-called viscous stirring which is proportional 
to the phase space average of A(e^) (see definitions of 
H\ and H 2 ). Depending on a er and &i T it can be either 
positive or negative. Second term represents the phe- 
nomenon of dynamical friction well known from galactic 
dynamics. In the limit m 3> m* its contribution is pro- 
portional to the mass of particle under consideration and 
to the surface mass density of field particles, but is inde- 
pendent of individual masses of field particles (cf . Binney 
& Tremaine 1987). This term is negative since it repre- 
sents the gravitational interaction of a moving body with 



the wake of field particles formed behind it as a result of 
gravitational focusing; thus H 2 < and K 2 < 0. Finally, 
the third term describes the increase of random velocities 
of a particular body at the expense of random motion of 
field planetesimals it interacts with. This effect is anal- 
ogous to the first order Fermi mechanism of the cosmic 
ray acceleration via scattering of energetic particles by 
randomly moving magnetic field inhomogeneities (Fermi 
1949). Note that in previous studies of planetesimal scat- 
tering it is the combination of the second and third terms 
on the r.h.s. of (JjjJ that is called the dynamical friction 
(Stewart & Wetherill 1988; Ida 1990). The combined ef- 
fect of these two terms is to drive planetesimal system to 
equipartition of random energy between planetesimals of 
different mass (which would be realized in the absence of 
viscous stirring). 

In this work we analyze planetesimal velocity evolution 
with the aid of equation ® . We call the first (positive) 
term on its r.h.s. gravitational stirring or heating (differ- 
ent from viscous stirring), while second (negative) term 
is called gravitational friction or cooling (different from 
dynamical friction 3 ). In this sense terms "stirring" and 
"friction" are only used to describe positive and negative 
contributions to the growth rate of planetesimals random 
motion. Use of equations © rather than © has two ob- 
vious advantages: (1) gravitational stirring and friction 
have different dependences on m* /(m + m*) which con- 
siderably simplifies analysis of velocity evolution equa- 
tions, and (2) stirring and friction terms have definite 
signs (unlike viscous stirring and dynamical friction used 
in previous studies of planetesimal dynamics). 

Gravitational scattering of planetesimals can proceed 
in two rather different velocity regimes, shear-dominated 
and dispersion-dominated. The former one is realized 
when a\ i + 0-* 2 <C (Rh/u) 2 , while the latter holds when 
(i^ + a* 2 ^> (Rh/o) 2 . Analytical arguments and nu- 
merical calculations demonstrate that in the dispersion- 
dominated regime o~ e and o~i are of the same order and 
evolve in a similar fashion (e.g. Ida 1990; Stewart & Ida 
2000). This happens because scattering in this regime 
has a three-dimensional character making different com- 
ponents of velocity ellipsoid comparable to each other. 
Thus, rough idea of the dynamical evolution in this case 
can be obtained by studying only one of the equations 
(J2J) . Moreover, the behavior of the scattering coefficients 
in the dispersion-dominated regime can be calculated an- 
alytically as a function of a er and di r in two-body ap- 
proximation (Stewart & Wetherill 1988; Tanaka & Ida 
1996), and one finds that 
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where Ai^Ci^ are positive functions of inclination to 
eccentricity ratio &i r /a er , and In A is a Coulomb loga- 
rithm (Binney & Tremaine 1987), which in our case can 
be represented as (Rafikov 2003a) 



A « aia ir (a 2 er + a 2 af r ), 



(7) 



3 I am grateful to Peter Goldreich and Re'em Sari for pointing 
this out to me. 



4 



with 0,1,2 being some constants. Explicit analytical ex- 
pressions for coefficients Ai. 2} Ci,2 have been derived by 
Stewart & Ida (2000). One can deduce a specific use- 
ful property of these functions by considering a single- 
mass planetesimal population: in dispersion-dominated 
regime the distribution of random energy between the 
vertical and horizontal motions tends to reach a quasi- 
equilibrium state (see e.g. Ida & Makino 1992). Then 
the ratio Oi r ja er is almost constant and both a er and a ir 
grow with time, meaning that [see © & (©] 

A x > 2A 2 , Ci > 2C 2 . (8) 

This property will be used later in JUfclCl 

In the shear-dominated regime o~ e and o~i evolve along 
different routes: eccentricity is excited much stronger 
than inclination because in this regime disk is geometri- 
cally thin and forcing in the plane of the disk is stronger 
than in perpendicular direction. Eccentricity evolution is 
also quite rapid in this case because of the vigorous scat- 
tering (relative horizontal velocity increases by ~ D,Rh 
in each synodic passage of two bodies initially separated 
by ~ Rh in semimajor axes). Simple reasoning con- 
firmed by numerical experiments suggests the following 
behavior of scattering coefficients in the shear-dominated 
regime: 
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where B\2,D\2 are some positive constants, which can 
be fixed using numerical orbit integrations, see Appendix 
El Simple qualitative derivation of © and can be 
found in Ida & Makino (1993) and Rafikov (2003c). 

For further convenience, we now reduce evolution equa- 
tions to a dimensionless form. We relate planetesimal 
masses to the smallest planetesimal mass 4 too by using 
dimensionless quantity x — m/rriQ. Instead of N(m) we 
introduce dimensionless mass distribution function f{x) 
such that 

v vj 

N(m) = -rf(x) and N{m)dm = -?-f(x)dx, (10) 
TOq mo 

where S p is the total mass surface density of planetesi- 
mals in the disk, which is a conserved quantity. We also 
rescale eccentricity and inclination dispersions: 
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(11) 



This is equivalent to rescaling all distances by Hill radius 
of smallest planetesimals. In this notation the boundary 
between the shear- and dispersion-dominated regimes is 
given by conditions 

s 2 + g *2 „ (a . + a .*)2/3 j fi 2 + s *2 ^ {x + a ,* ) 2/3_ (J2) 

We also introduce dimensionless time r by 

-2/3 



t 



where to 
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4 Protoplanetary disks should contain planetesimals of different 
sizes but we assume that bodies lighter than mo are not important 
for the disk evolution. In the case considered here the choice of mo 
is dictated by the mass at which distribution would depart from 
the given power law form towards shallower scaling with mass (see 
also §4.1.1). 



is a timescale for the stirring of planetesimal disk com- 
posed of bodies of mass too only on the boundary be- 
tween the shear- and dispersion-dominated regimes. Nu- 
merically, for M c = Mq one finds that 

^ lOgcm- 2 ^ / m \ 1/3 
t w 16 yr a A u ^ 77^7- , (14) 



where we assumed that Tip (a) = T, p oa A fj (Hayashi 
1981). Dynamical timescale is rather short when s, s z ~ 
1, but it grows very rapidly as epicyclic velocities of plan- 
etesimals increase (t ~ tos 4 in the dispersion-dominated 
regime). 

Random velocities of planetesimals of mass x evolve 
due to the interaction with all other bodies spanning the 
whole mass spectrum. Using our dimensionless notation 
we can rewrite equations (0, ©, & © in the following 
general form: 
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where s = s(x), s* 
shear-dominated regime [s 2 

H\ » Ci, H 2 ~ — C2 



x* sj + si 

s(x*) and similarly for s z , s*; in the 
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and in the dispersion-dominated regime 
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(18) 

s + s 

We explore this system in 21 using two approaches: 
asymptotic analysis utilizing analytical methods, and di- 
rect numerical calculation of velocity evolution. 

3. MASS SPECTRUM. 

Starting from (|15l) - (|18fl one would like to obtain the 
behavior of s and s z as functions of x and r, given some 
planetesimal mass spectrum f{x) as an input. In gen- 
eral, to do this one has to evolve equations i|16|) - (|18fl nu ~ 
merically. However, it is usually true that the planetes- 
imal size distribution spans many orders of magnitude 
in mass; thus one would expect that some general pre- 
dictions can be made analytically about the asymptotic 
properties of planetesimal velocity spectrum. 

In this study we assume that planetesimal mass distri- 
bution has a "self-similar" form; specifically we take 



f(x,r) = ip{r)ip 



x c (r) 



(19) 
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Fig. 1. — Schematic representation of the shape of planetesimal 
mass spectrum and its evolution in time. Mass distribution has a 
power law form exponentially cut off at x c (t) with a tail of runaway 
bodies for x x c (t). 

where x c (t) 3> 1 is some fiducial planetesimal mass 
which steadily grows in time as a result of coagulation, 
V>(t) is a temporal modulation, and function ip repre- 
sents a self-similar shape of the mass spectrum. We also 
assume in this study that asymptotically 

<p(v)~l v ~ a l ^ y *b (20) 

^ yy > \exp(-y), y ' 

where a > 0. Such scaling behavior is often found in 
coagulation simulations. We will see in 21 that for clari- 
fying the asymptotic properties of planetesimal velocity 
spectrum it is enough to know only the asymptotic be- 
havior of (p(y) and not its exact shape. Moreover, in |B] 
we demonstrate how our results for power-law size dis- 
tributions can be generalized for other mass spectra. 

Normalization iP(t) is not an independent function. It 
is related to x c (t) because of the conservation of the total 
(dimensionless) planetesimal surface density 

oo 

Mi = J xf(x)dx. (21) 
l 

Taking this constraint into account one finds that 

In the case a > 2 mass spectrum behaves as f(x,r) — 
x~ a for x <C x c (t); only the position of the high mass 
cutoff x c (t) shifts towards higher and higher masses with 
time. 

We define M v (t) to be the ^-th order moment of the 
mass distribution: 

oo 

M v {r) = j{x*Yf{x\r)dx\ (23) 
l 



Integral in (|23() is dominated by the upper end of the 
power-law part of the mass spectrum (i.e. by x* ~ x c ) if 
v > a — 1; in this case we can write [using (|19|) ] 

M v {t) =M v [x c {r)Y +1 i){r) 

o 

where M v — reduced moment of order v — is a time- 
independent constant for a given mass distribution. 

One of the interesting features exhibited by self- 
consistent coagulation simulations is the development of 
the tail of high mass bodies beyond the cutoff of the 
bulk distribution of planetesimals (Wetherill & Stewart 
1989; Inaba et al. 2001). This is interpreted as the man- 
ifestation of the runaway phenomenon taking place in a 
coagulating system. To explore the possibility of the run- 
away scenario we add to the mass spectrum (12011 a tail of 
high mass planetesimals extending beyond the exponen- 
tial cutoff x c (see Appendix lA"|l. In doing this we always 
make sure that runaway tail contains a negligible part 
of the system's mass and does not affect its dynamical 
state (i.e. stirring and friction caused by these high mass 
planetesimals are small) which is true if biggest bodies 
are not too massive (Rafikov 2003c). Under these as- 
sumptions the explicit functional form of the runaway 
tail is completely unimportant. Throughout this study, 
we use the term "planetesimals" for bodies belonging to 
the power-law part of the spectrum (x < x c ), and "mas- 
sive" or "runaway" bodies to denote the constituents of 
the tail (x ^> x c ). 

We do not restrict the variety of possible functional 
dependences of x c on time r. It will turn out that all our 
results can be expressed as some functions of x c which 
leads to the time-dependence of planetesimal velocities 
in a very general form. In some cases it will be impor- 
tant that x c (t) does not grow too fast, which, however, 
we believe is a fairly weak constraint (see discussion in 
§4.1.2). 

4. VELOCITY SCALING LAWS. 

To facilitate our treatment of planetesimal velocities 
we introduce a set of simplifications into our considera- 
tion: 

• We split mass spectrum into several regions, such 
that in each of them planetesimal interactions 
can be considered as dominated by a single dy- 
namical regime (shear-dominated or dispersion- 
dominated). Transitions between such regions are 
not considered but in principle might be treated by 
interpolation. 

• When studying the dispersion-dominated regime 
we neglect the difference between the vertical and 
horizontal velocity dispersions, and treat them by 
a single equation. We also set Coulomb logarithm 
equal to constant, because of its rather weak de- 
pendence on s, s z or x. 

The validity and impact of these assumptions on the ve- 
locity spectrum are further checked using numerical tech- 
niques. Since we are mostly interested in qualitative be- 
havior of solutions, we do not expect numerical constants 



() 



to be correctly reproduced by our analysis; they can be 
trivially fixed using numerical calculations. 

Mass distributions are parametrized by the value of 
power law exponent a [see (|20fl ]. We split our investi- 
gation according to the value of a into 3 cases: shal- 
low mass spectrum, a < 2, intermediate mass spectrum, 
2 < a < 3, and steep spectrum, a > 3 (this regime splits 
into two more important subcases, see §4.3). Note that 
for the sake of avoiding additional complications we do 
not consider borderline cases a = 2,3; they can be eas- 
ily studied in the framework of our approach if the need 
arises. 

In our analytical work planetesimals are started with 
large enough s and s z so that they interact with each 
other in the dispersion-dominated regime. We then 
also assume that planetesimals with masses x < x c (r) 
(containing most of the mass) stay in the dispersion- 
dominated regime w.r.t. each other also at later time, 
i.e. 

s(x,t), s z {x,t) > x 1/3 for x<x c (t). (25) 
This is a reasonable assumption for all mass spectra at 
the beginning of evolution, although for steep size distri- 
butions it may break down when maximum planetesimal 
mass becomes very large (see §4.3). 

4.1. Shallow mass spectrum. 

Planetesimal mass distributions shallower than m~ 2 
result from coagulation of high-velocity planetesimals, 
when gravitational focusing is unimportant and collision 
rate is determined by geometrical cross-section of collid- 
ing bodies (Wetherill & Stewart 1993; Kenyon & Luu 
1998). It is worth remembering however that in highly 
dynamically excited disks (1) energy dissipation in in- 
elastic collisions must be important, see SJSJ an d (2) plan- 
etesimal fragmentation cannot be ignored. Despite that, 
the case of shallow mass spectrum is interesting because 
it facilitates understanding of the disks with other plan- 
etesimal size distributions. 



4.1.1. Velocities of planetesimals. 
We start by considering planetesimals, x < 



x c (t) 



part of the size distribution containing most of the 
mass. Regarding them as interacting in the dispersion- 
dominated regime [i.e. the condition (|25|l is fulfilled] we 
may write using (flt)|l & (fl%|l that 



0s* 
~d7 



dx*f{x*) 



-Ai - 2- 



x*(x + X*) 



o*2 



„*2 



A 2 



(26) 

_x + X* s* + S** 

for x <C x c {t). Since we use a single equation to describe 
the evolution of both s and s z , the values of constants 
are not well defined but this is not important for 
deriving general properties of the velocity spectrum. 

We can represent (|26fl asymptotically in the following 
form: 



ds^ 
dr 



,.*2 



A x \ dx*- 



/(**) 



„*2 



dx*~- 



/(**) 



At 



2A 2 s 2 x / dx* 



x*f(x*) 
(s 2 + s* 2 ) 2 



(27) 



S,S 




shear 



Fig. 2. — Sketch of the behavior of the scaled planetesimal eccen- 
tricity and inclination s = a £ a/ru and s z = oia/rn as functions 
of their dimensionless mass x = m/m,Q, predicted by asymptotic 
analysis of §4.1. Solid line displays s(x) while dashed line is for 
s z (x) (for x < Xghear their behaviors are the same). Dotted line 
delineates shear- and dispersion-dominated regimes of planetesimal 
interaction with bodies of comparable mass. The rest of notation 
and theoretical predictions for the temporal behavior of s(x) and 
s z (x) can be found in the text. 



The first two terms on the r.h.s. represent correspond- 
ingly the gravitational stirring and friction by bodies 
smaller than x. The last term is responsible for the com- 
bined effect of stirring and friction produced by bodies 
bigger than x. From this equation it is easy to sec that 
s(t) independent of x is a legitimate solution of (|27ll in 
the case a < 2. Indeed, assuming s to be independent of 
x we find that all integrals over x* in (|27J) are dominated 
by their upper limits when mass spectrum is shallow. 
Using (|19fl - H24fl we can estimate that 



ds 2 ^ ^ x c ix 



dr 



2s 2 



A 2 ) 



M 2 (r) 
2s 2 ' 



(28) 



Note that in deriving this expression we have extended 
the integration range of the last term in l|27|l from 1 to 
infinity (not from ~ x). This is legitimate because this 
integral is dominated by the contribution coming from 
~ x c 3> x and, thus, adding interval from 1 to ~ x to 
integration range would introduce only a subdominant 
contribution to the final answer. Clearly, two first terms 
in l|28|l are negligible compared to the third one for x -C 
x c , and third term is positive [see JSJ] and independent 
of x, meaning that s(r) is also independent of mass, in 
accord with our assumption. One can then easily find 
that 5 [dropping constant factors A% t 2 but bearing in mind 

5 Throughout the paper we imply by so the velocity dispersion of 
sma llest plan etes imals and so is different for different mass spectra, 
see 141)1 and 1531 below. 
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condition 



x. x„ 



s(x,t) w s 



M 2 (t')cIt' 



1/4 



M 2 / x c (r')dr' 



1/4 



a < 2, x < x c (29) 



[if s(x,t) 3> s(x, 0)]. The second equality holds only 
for a < 2, according to (22) & (ED- Thus, in the case 
of mass spectrum shallower than x~ 2 random velocities 
are independent of planetesimal masses and uniformly 
grow with time. Both gravitational stirring and friction 
are dominated by biggest planetesimals (x ~ x c ), with 
stirring being more important; friction by small bodies is 
completely negligible and energy equipartition between 
planetesimals of different mass is not reached. For this 
reason, planetesimal velocity expressed in physical units 
using Hll|) , 1|13|) , & (|29l) is independent of the choice of 
m . A schematic representation of velocity spectrum for 
the shallow mass distribution is displayed in Figure [21 

4.1.2. Velocities of runaway bodies. 

Now we turn our attention to the runaway bodies, x > 
x c . By assumption, this tail contains so little mass, that 
it cannot affect velocities not only of planetesimals but 
also of runaway bodies themselves. 

Because of our assumption of the dispersion-dominated 
interaction between all planetesimals with x < x c it 
is natural to expect that at least those runaway bod- 
ies which are not too heavy still experience dispersion- 
dominated scattering by planetesimals with masses < x c . 
Assuming that the velocity dispersion of these runaway 
bodies is much smaller than that of the bulk of planetes- 
imals — natural consequence of the tendency to equipar- 
tition of random epicyclic energy between bodies of dif- 
ferent mass — we can rewrite <|26l) as 



9s 2 „ 



dx*- 



„*2 



- 2A 2 s 2 x 



Although the integration range of all terms in the r.h.s. 
is extended to infinity, the exponential cutoff of f(x) at 
~ x c effectively restricts the integration range to be from 
1 to ~ x c . 

One can easily see that s[x, r) = const(x) cannot be a 
solution of (|30(l because dynamical friction is too strong. 
Indeed, if we were to assume the opposite, we would find 
that s(x, t) is still given by (|29|l . However, direct substi- 
tution of (|29|l into (|3(J[1 shows that the gravitational fric- 
tion term exceeds other contributions for x > x c . This 
might tempt one to suggest that heating (first) term in 
r.h.s. (|3U[) should be neglected compared to the friction 
(second) term. In this case one would find that s(x,t) 
exponentially decays in time and very soon heating term 
becomes the dominant one, contrary to initial assump- 
tion. These results suggest the only remaining possibil- 
ity — that heating of massive bodies by planetesimals 
almost exactly balances friction by planetesimals, and 
minute difference between them is accounted for by the 
time dependence of s(x, r) embodied in the l.h.s. of (|30|) . 




2 4 



10 12 14 



Log(x) 



Fig. 3. — Numerical results in the case of a = 1.5 for (a) ec- 
centricity dispersion (dimensionless) s and (b) inclination disper- 
sion s z . In the bottom panels the velocity spectrum is plotted as a 
function of dimensionless mass x, while in the top panels the power 
law index r\ = din s/ din x of the spectrum is displayed. Different 
curves correspond to increasing values of maximum planetesimal 
mass x c (and, consequently, time) indicated by arrows of corre- 
sponding line type on the upper panels. Dotted line in the lower 
panels is s = x 1 / 3 (boundary between different velocity regimes). 
Dotted curve in the upper panels is the theoretical prediction for 
the run of the spectral index rj corresponding to the largest x c dis- 
played (for which numerical result is always shown by thick solid 
line). These results are to be compared with theoretical predictions 
shown in Figure 121 



Given that this assumption is correct we immediately 
find that 



s(x,t) 



dx 



^* 2 f(x*) 

„*2 



1/2 



dx 



s *4 



-1/2 



(31) 



It is easy to see that scaling s oc x 1//2 follows di- 



rectly from assuming that (1) the interaction is in the 
dispersion-dominated regime, and (2) heating of runaway 
bodies by planctcsimals is balanced by friction also due 
to planetesimals. This specific result is completely inde- 
pendent of the mass or velocity spectra of small bodies. 
We will highlight this again when we obtain similar re- 
sults for intermediate and steep mass spectra in §4.2.2 
and §C.2. For a specific case a < 2 considered in this 
section s* — sq(t) and one finds that 

"M 2 (r)l 1/2 



s(x,t) 



So(t) 




< x < 



%shear 



(32) 



[x shear is defined below in equation JHHJ] . At x ~ x c this 
formula is consistent with so — extrapolation of small 
mass result (|29H up to x ~ x c . 

Time dependence enters (|32|) eventually through x c (r). 
Thus, it is the behavior of x c (t) that determines ds 2 /dr 
in The timescale on which x c varies is usually much 
longer than the dynamical relaxation time of the system, 
thus one would expect l.h.s. of l|30(l to be small in ac- 
cordance with the assumptions which led to (|31|l . The 
requirement of the negligible effect of time derivatives in 
situations like the one described here is the only con- 
straint which must be imposed on the possible behavior 
of x c (t). 

As we consider even bigger x we find that at some mass 
runaway bodies start interacting with bodies of similar 
mass in the shear-dominated regime; this happens when 

x 1 / 3 ~ s (x c /x) 1 / 2 , or x ~ x c (s /xc /3 ) 6/5 . This, how- 
ever does not affect the dynamics of these bodies, be- 
cause by our assumption runaway tail contains too little 
mass to dynamically affect its own constituents. It is 
only after the runaway bodies start interacting in the 
shear-dominated regime with planetesimals (i.e. bod- 
ies lighter than x c ), that the velocity evolution of the 
tail gets affected. In our case this clearly happens when 
x 1 / 3 ~ sq(t). Thus, runaway bodies with masses 

X ^ X s /j ear = Sq (33) 

interact with all planetesimals in the shear-dominated 
regime. 

As mentioned in |J21 in the shear-dominated case ec- 
centricities and inclinations may no longer evolve along 
similar routes and we have to consider them separately. 
Using lO & (O we write 



dr 
Ci 



dx*f(x*)x*(x + x*) 1 / 3 



l 

x* 



x + x* 

oo 



2C 2 



(x + x*) 2 / 3 
ds 2 f s 2 + s* 2 



(34) 



)l/3 



' x + x* 



-2D 2 - 



(35) 



Assuming that s, s z -C s* ~ s* (because of the gravi- 
tational friction) and x 3> x* (both heating and cooling 
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Fig. 4. — Comparison of numerical (solid curves) and analytical 
(dotted curves) behaviors of sq as a function of dimensionless time 
r. We display cases: (a) a = 1.5, (b) a = 2.3, (c) a = 2.7, (d) 
a = 3.2 and a = 4.3. 



are dominated by planetesimals, not runaway bodies), 
and neglecting l.h.s. of both equations (similar to the 
previously discussed situation for x c < x < x s } lear ), we 
obtain the following result: 



s(x,r) » x- 1 / 6 



s z (x,t) w — 
VX 



M 2 (r) 



-,1/2 



Mi 



oo 

/ dx*x* 2 s* 2 fix*) 

Mi J V ' 



1/2 



(36) 



Finally, setting s* sq(t) as appropriate for a < 2 we 
find that 



s(x, t) w x 



-1/6 



M2 

Mi 



1/2 



S ,(x,r)« S0 (-j — 



1/2 



x ^ x s h ear . (37) 



These expressions demonstrate that transition to the 
shear-dominated regime is accompanied by a change of 
the power-law index of eccentricity scaling with mass, 
while the power-law slope of inclination scaling stays 
the same. This happens because in the case of inclina- 
tion evolution heating term is functionally similar to the 
friction term (exactly like in the dispersion-dominated 
case), which is not the case for eccentricity evolution, see 
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<|17|) . As a result, for x < x shear one finds that s 3> s z 
— planetesimal velocity ellipsoid flattens considerably in 
the vertical direction, which is very unlike the dispersion- 
dominated case (see Figure 0). This conclusion is very 
general and should be found whenever runaway bodies 
interact with all planetesimals in the shear-dominated 
regime. 

4.1.3. Comparison with numerical results. 

To check the prediction of our asymptotic analysis 
we have performed numerical calculations of planetesi- 
mal velocity evolution for a distribution of planetesimal 
masses which was artificially evolved in time according 
to prescriptions outlined in Details of our numerical 
procedure can be found in Appendix 1X1 

In Figure we show the evolution of planetesimal ec- 
centricity and inclination dispersions for a shallow mass 
spectrum with a = 1.5. For each curve displaying s or 
s z run with x at a specific moment of time we also com- 
puted power law spectral index 77 = alhis/dliix. We 
then compare power law indices r\ and r\ z with analyt- 
ical predictions given by equations {Sjj), & (|37|) . 
Dotted line in the upper panels of Figure represents 
analytical predictions for 77 and r\ z at the moment when 
the snapshot of velocity spectrum displayed by the solid 
line (corresponding to the highest x c shown) was taken. 
One can see that agreement between two approaches is 
quite good. Power law slopes of both eccentricity 77 and 
inclination r/ z exhibit a well defined plateau at masses 
considerably smaller than x c which is predicted by (1291) . 
Above upper mass cutoff 77 and r\ z plunge towards —1/2 
which is in agreement with l|32[l . although they do not fol- 
low this value very closely. The reason is most likely the 
lack of dynamical range — transitions between different 
regimes typically occupy substantial intervals in mass. 
Finally, at x above x s } lear , where runaway bodies inter- 
act with all planetesimals in shear-dominated regime, 77 
goes up to —1/6 while r\ z stays at a value of —1/2, ex- 
actly like equation (|37[1 predicts. Glitches in the curves 
of 77 and r\ z in the vicinity of x s h ea r are the artifacts of in- 
terpolation of scattering coefficients between shear- and 
dispersion-dominated regimes. 

In Figure^ we display the time evolution of so — ec- 
centricity dispersion of smallest planetesimals versus the 
prediction of our analysis given by equation 1(29(1 . One 
can see that the agreement between two approaches is 
excellent — analytical and numerical results closely fol- 
low each other as time increases. Vertical offcet between 
them amounts to a constant factor by which numerically 
determined so(t) differs from analytical sq(t). This, of 
course, is expected because we completely ignored con- 
stant coefficients in our asymptotic analysis. Numerical 
results shown in Figure 0]} allow us to fix constant coef- 
ficient in (|29|l . and we do this when we calculate x s h ear 
in Figure using equation 1)33(1 . Note that the rapid 
growth of so with time should certainly increase veloc- 
ities of planetesimals beyond their escape speeds (if it 
were not the case initially). Dynamical evolution in this 
case is discussed in more details in 

An apparent break in the behavior of so(t) at r ~ 10 5 
can be seen in Figure 0^. It occurs because this is the 
time at which mass scale x c starts to grow. The pre- 
scription for x c (t) which we use (described in Appendix 
lAl is such that x c is almost constant for r < 10 5 ; as a 
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Fig. 5. — Same as Figurel2*lbut for the case of intermediate mass 
spectrum. Specific case of2<«<5/2is displayed. Notation is 
explained in the text. Analogous plot for 5/2 < a < 3 would 
exhibit different behavior of s z in the high-mass tail and can be 
easily constructed using the results of §4.2.2. 



result, so <x t 1 / 4 according to equation J2HJ). At r > 10 5 
scale factor starts to increase rapidly and this causes so 
to grow faster. To conclude, the results of this section 
show that the overall agreement between numerical cal- 
culations and analytical theory of planetesimal velocity 
evolution is rather good. 

4.2. Intermediate mass spectrum. 

Self-consistent coagulation calculations (Kenyon & 
Luu 1998) and N-body simulations of planet formation 
(Kokubo & Ida 1996, 2000) often yield intermediate mass 
distributions of planetesimals (2 < a < 3). Kokubo & 
Ida (1996) have found, for example a w 2.5 ± 0.4 in 
their N-body calculation of collisional evolution of several 
thousand massive (m = 10 23 g) planetesimals. Such out- 
come seems to be ubiquitous for planetesimals which ag- 
glomerate under the action of strong gravitational focus- 
ing (planetesimal velocities are well below the escape ve- 
locities from their surfaces) in the dispersion-dominated 
regime, which makes the case of intermediate mass spec- 
tra very important. 

4.2.1. Velocities of planetesimals. 

As in §4.1 we start our consideration of mass spectra 
with 2 < a < 3 by concentrating on planetesimals with 
masses < x c (t), assuming that all these bodies interact 
in the dispersion-dominated regime (see f^J. Then equa- 
tions Q26J) and 1)27(1 continue to hold. 

It is natural to start by assuming again that s is in- 
dependent of x. Performing analysis similar to that of 
§4.1.2 we find that stirring is still dominated by large 
planetesimals (~ x c ), but the biggest contribution to the 
gravitational friction is produced by smallest planetes- 
imals (with to ~ toq, or i ~ 1 in our dimensionless 
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notation) : 
dr 



„3-a 



2s 2 



2A 2 — Mi 



(Ax - A 2 



M 2 (t) 



(38) 



2s 2 

Clearly, the first term on the r.h.s. which is responsi- 
ble for the stirring by small planetesimals is negligible 
compared to the last one, representing stirring by plan- 
etesimals bigger than x [see Q24jl]. Also, by definition, 
Mi is just a dimensionless total surface density of plan- 
etesimal disk, and thus has to be constant in time. Then 
our assumption of s independent of x is self-consistent 
and s is given by the first equality in (|29|l only if friction 
by small planetesimals [ second term in (|38(l ] is small 
compared to stirring by big ones (third term), which is 
true provided that 



x<x k (r) = — 



M2 2 . 
'M 1 Xc 



(39) 



The last equality follows from the fact that M 2 = 
M2X^T a when 2 < a < 3. Thus, for small planetesimals 
(x < %k) friction by even smaller ones is again unim- 
portant; their velocities are excited by big planetesimals 
(x r~j x c ), and as a result [cf. I|29|) ] 

1/4 



s(x,t) ps s (r) 



M 2 / [x c { T ')f- a dT' 



2 < a < 3, x < x k . (40) 
t.u. gravitational friction by small 



= const(x), 

However, for x £ Xk, 
planetesimals is no longer negligible compared to the ve- 
locity excitation by big ones. One would then expect 
friction to lower the velocities of big bodies below those 
of small planetesimals (because of the tendency to en- 
ergy equipartition associated with friction term). Start- 
ing with (|27|) and using 1401) we obtain the following equa- 
tion for x > Xk'- 



S" = l| / dx^x^fix^ + A, J 



dx 



*(z*) 2 /(z*) 



,*2 



s 4 
• s 



2A 2 — / dx*x*f(x*)-2A 2 s 2 x / dx* 
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Fig. 6. — The same as Figure[3]but for the case of intermediate 
mass spectrum with a = 2.3 (cf. Fig. |SJ. 



A x - 2A 2 



dx*{x*) 2 f{x*) 



(41) 



In this equation different terms represent stirring or fric- 
tion by different parts of planetesimal mass spectrum, 
below and above Xk- We now take an educated guess 
that the solution for s can be obtained by equating third 
and fifth terms in the r.h.s. of (|41|l . Physically this means 
that velocity spectrum results from the balance of veloc- 
ity stirring by biggest planetesimals (masses ~ x c x, 
similar to the case of shallow mass spectrum) and friction 
which is mainly contributed by smallest planetesimals. 
One then finds that [see Q39[l] 



s(x,r) 



»o(t) 



so(t) 

Xkij) 



M 2 (t) 
Mi 

-,1/4 



1/4 



.-1/4 



Xk 



<T< 



(42) 



Substituting this solution into (|41|l we easily verify that 
our neglect of all other terms in the r.h.s. of (|41|) is 
indeed justifiable. The l.h.s. of (|4*T)> can be neglected for 
x c on the basis of arguments identical to those 



x k 



<x< 



advanced in §4.1.2. Thus, s given by (|42|l is indeed a 
legitimate solution in this mass range. 

Equations l|40|l & (|42|l highlight an interesting feature 
of the case 2 < a < 3: although the mass spectrum 
has a power law form with a continuous slope all the 
way up to x c , velocity spectrum exhibits a knee at some 
intermediate mass Xk (see Figure [EJ . The resemblance 
of the velocity spectrum for x < Xk to the spectrum for 
a < 2 is due to the fact that in both cases stirring is done 
by biggest planetesimals and friction has a small effect. 
However, for x > Xk effect of friction becomes important, 
and now it is dominated by the part of mass spectrum 
(smallest planetesimals) different from that responsible 
for stirring (biggest planetesimals). This causes velocity 
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spectrum to change its slope from to —1/4 in this mass 
range (see Figure EJ. 

4.2.2. Velocities of runaway bodies. 

Similar to the case of shallow mass spectrum, some 
(lightest) runaway bodies interact with the bulk of plan- 
etesimals in the dispersion-dominated regime. In this 
case all the considerations of §4.1.2 hold, i.e. evolution 
equation (|3U|> stays unchanged and after the same line of 
arguments we arrive at the expression l|31|l for the plan- 
etesimal velocity dispersion. Integrals entering (|31|) must 
be reevaluated anew using equations (|40|l & (|42|l appro- 
priate for 2 < a < 3. Performing this procedure and 
substituting resulting expressions into (|31fl we find that 



S[X,T 



so 
Jx 



M 5/2 (t) 



1/2 



so 



(?) 



1/2 / Mr 




This velocity distribution is valid only for runaway bodies 
lighter than 



Xh 



[s(x c )Y 



3/4 



(44) 



because only such bodies interact with all small mass 
planetesimals (x < x c ) in the dispersion-dominated 
regime. In agreement with what we found earlier in 
§4.1.2, velocity dispersion decreases with mass as x^ 1 !" 1 . 
At x <~ x c this expression matches the low mass result 

Runaway bodies heavier than X] x but still lighter 
than x s hear = Sq are in a mixed state: they in- 
teract with the most massive planetesimals in the 
shear-dominated regime, and with lighter ones in the 
dispersion-dominated regime. This, of course, is a di- 
rect consequence of the fact that the velocity spectrum 
of planetesimals has a knee at Xk, see §4.2.1. Introducing 



(45) 



we find that runaway bodies with masses x between 
Xh and x s hear interact with planetesimals less massive 
than x s (x) in the dispersion-dominated regime, and with 
planetesimals heavier than x s {x) in the shear-dominated 
regime. Self-consistent analysis of planetesimal velocity 
spectrum for Xh < x < x s hear is described in detail in 
Appendix [5] and only final results are shown here. It 
turns out that eccentricity scales as 
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(46) 



while inclination behaves in a different manner: 
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Fig. 7. — The same as Figure[3]but for the case of intermediate 
mass spectrum with a = 2.7 (cf. Fig. |SJ. 



s z (x,t) PS s 

„ „(4a-17)/6 



( Xs (x)y/ 2 - a 

x 1 J 2 xM l 



1/2 



a>5/2, (48) 

In the case of a = 5/2 a more general expression for s z 
following directly from l|B4|l should be used. 

Finally, bodies heavier than x s hear interact with all 
planetesimals in the shear-dominated regime. Expres- 
sions l)36|l adequately describe velocity behavior in this 
case. Manipulating them with the aid of equations l|4l)|l 
and (|l2"jl we find that 



s{x,t) 



1/2 
,1/6' 



^ & shear • 



(49) 



while 



a < 5/2, (47) s z (x, r) « s (^) ^ (^) ^ 



3/2 



1/2 
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s z (x,t) sa s 



x) 



(3-a)/2 



a < 5/2, (50) 
> 5/2. (51) 



Scalings of velocity dispersions with mass x in (|49|) - (|51|l 
are the same as in the case of shallow mass spectrum, but 
the time dependences are different. Velocity behavior in 
the case of intermediate mass spectrum is displayed in 
Figure El 

4.2.3. Comparison with numerical results. 

We calculated numerically the evolution of s and s z for 
two intermediate mass spectra: a — 2.3 and a = 2.7. In 
Figure 0d,c we compare numerically calculated sq(t) - 
eccentricity dispersion of smallest planetesimals — with 
predictions of our asymptotic theory. As in the case of 
shallow mass spectrum discussed in §4.1.3, agreement be- 
tween the two approaches is very nice. This allows us to 
fix coefficient in l|40|l , absence of which in our asymptotic 
analysis causes a constant offset between numerical and 
analytical curves in Figure 0Jd,c. Results for the velocity 
dependencies on mass at different moments of time and 
their comparison with analytical theory of intermediate 
mass spectra outlined above are shown in Figures El and 

In the case a = 2.3 one can clearly see the appearance 
of all features we described in §4.2.1 and §4.2.2. Indeed, 
in Figure El one can observe both the zero-slope part of 
the spectrum at x < x k and the trend of reaching the 
slope —1/4 for x k < x < x c ; the last feature is not 
fully developed when we stop our calculation but it is 
almost certain that rj and r\ z would reach slope —1/4 
given enough time and mass range. Calculation has to 
be stopped when x reaches ~ 10 8 because later on biggest 
planetesimals start interacting with bodies of similar size 
in the shear-dominated regime. Such a possibility was 
not explicitly considered in the present study (but it can 
be easily dealt with using analytical apparatus developed 
in this work). 

Results for biggest runaway bodies (x > x shear) are 
in nice concordance with analytical predictions given by 



aquations (|51|l . In the range x c < x < s shear velocities of 
the runaway bodies clearly have not yet converged to a 
steady state but are evolving in the right direction. The 
range of masses where velocity spectrum should exhibit 
a slope —l/2[x c <x<Xh, see (|43|) ] is so narrow 6 for the 
last (solid) curve (because planetesimals with mass ~ x c 
are already very close to shear-dominated regime) that 
we do not see this regime realized. At the same time, it 
is conceivable that runaway bodies which interact with 



< X < 



3 she 



see 



planetesimals in the mixed regime (xh 
Appendix will finally reach their asymptotic velocity 
state with n = —5/6 and n z — —7/6 [see equation 
although it will take very long time for them to get there. 

Comparison for the case a = 2.7 (Figure 01 is very 
similar. All major features corresponding to the interme- 
diate mass spectrum can be traced in this case as well. 
Unfortunately, the comparison with analytical results is 
complicated by the fact that in this case we have to stop 

6 Position of in Figure[S]is determined by equation 1441 which 
assumes that between x k and x c velocity behaves oc x~ 1//4 . In 
reality, as we discussed abo ve, velocity slope in that mass region is 
somewhat steeper causing 1441 to overestimate x^. 
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Fig. 8. — (a) Same as Figurel2*lbut for the Case 1 of steep mass 
spec trum (3 < a < 4). Note that 1/4 < /3 < 1/2 [see equation 
1551 1 . (b) Same as Figure |5] but for steep mass spectrum with 
a > 4. 



numerical calculation even earlier than in a = 2.3 case 
to avoid bringing most massive planetesimals into the 
shear-dominated regime (for this reason mass scale in 
Figure does not go as far as in Figure El • As a re- 
sult, even a velocity plateau at small masses does not 
have time to fully develop, not speaking of runaway bod- 
ies in "mixed" interaction state. Despite this, the overall 
agreement between the numerical results and asymptotic 
theory of velocity evolution of intermediate mass distri- 
butions is rather good, especially if we make allowance 
for the short duration of our calculation and small mass 
range of planetesimals. 

4.3. Steep mass spectrum. 
Planetesimal spectra decreasing with mass steeper 



than 



,-3 



are usually not found in calculations of plan- 



etesimal agglomeration. Still, this case is quite interest- 
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Fig. 9. — The same as FigurelSlbut for the case of steep mass 
spectrum with a = 3.2 (cf. Fig. 18b). 



ing and we briefly discuss it here. 

As demonstrated in §4.1 & §4.2 stirring of planetesi- 
mals is determined solely by biggest bodies with masses 
~ x c when a < 3. This is no longer true when the mass 
spectrum becomes steeper than x~ 3 . Indeed, let us con- 
sider the evolution of sq — velocity dispersion of smallest 
planetesimals having dimensionless mass x = 1. Using 
equation (|26[) we find that 



34 
dr 



1 



dX I{X ] 1 + s^/sl 



-A 1 



1 + s^/sl 



(52) 



Assuming that s*/s — ► when x* fx — > oo as a re- 
sult of gravitational friction we find that integral in 152f) 
is dominated by its lower integration limit for a > 3. 
This means that the velocity of smallest planetesimals is 



now mediated by smallest planetesimals themselves, un- 
like the case of a < 3. We also find that the temporal 
behavior of sq is given simply by 



s (r) 



,1/4 



a > 3. 



(53) 



We now turn our attention to studying planetesimal 
velocities in the mass range 1 < x < x c . We make an 
a priori assumption that planetesimal velocity spectrum 
has a form of a simple power law (with constant power 
law index) 

s(x, t) « s (t)^ _/3 , x < x c , (54) 

and verify later if it is correct. In Appendix IH1 we show 
that whenever 3 < a < 4, both heating and friction of 
planetesimals of a particular mass x are determined by 
bodies of similar mass, x* ~ x. This is somewhat unusual 
in view of our previous results when only the extrema of 
power-law mass spectrum were important (see §4.1 and 
§4.2. In this case power law ind ex 8 can only be found 
by numerically solving equation (|C2(I : there is no simple 
expression for 8 but its value has to satisfy condition 



4 2 



(55) 



For even steeper mass distributions, a > 4, we demon- 
strate in Appendix IU1 that smallest planetesimals domi- 
nate both stirring and friction not only of themselves, 
but also of bigger bodies. This occurs because mass 
spectrum is very steep and high mass planetesimals con- 
tain too little mass to be of any dynamical importance. 
This is very similar to the velocity evolution of the run- 
away tail in the dispersion-dominated regime described in 
§4.1.2, and, not surprisingly, we find that 8 = 1/2 when 
a > 4. Thus, a solution leading to a complete energy 
equipartition between planetesimals (s cx x~ 1 / 2 ) is only 
possible for very steep mass spectra. Velocity behavior 
of runaway bodies in the case of steep mass spectrum is 
described in detail in Appendix [U] Theoretical spectra 
of planetesimal velocities are schematically illustrated in 
Figures (for 3 < a < 4) and|Hb (for a > 4). 

We verify the accuracy of these predictions by numer- 
ically calculating the velocity evolution of planetesimal 
disks with a = 3.2 and a = 4.3. In the first case, shown 
in Figure|^l the slope 8 of planetesimal velocity spectrum 
for x < x c is evidently shallower than —1/2 but steeper 
than —(a — 2)/4 — —0.3, in agreement with constraint 
(|55|) . Analytical prediction for the behavior of rj(x) in 
top panels of Figure is computed using the numeri- 
cally determined value of 8 and equations l|C4|l . (|C5(I . 
and i|C6fl . It agrees with numerical result (solid curves 
corresponding to the largest x c displayed) rather well, 
especially for x > x shear- 

In the case a = 4.3 (Figure I10[l the slope of plan- 
etesimal velocity scaling with mass is essentially indis- 
tinguishable from —1/2 for x < x c , in accordance with 
our asymptotic prediction. Velocity evolution of runaway 
tail also agrees with the discussion in Appendix C.2 quite 
well. In both a = 3.2 and a — 4.3 cases the final velocity 
curves displayed are at the limit where we can still apply 



our theory: s(x c 



,V3 



; this washes out some details of 



velocity spectra of runaway tails predicted in Appendix 
C.2. Differences of r] and r] z from their predicted values 
at smallest planetesimal masses are caused by boundary 
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Log(x) 

Fig. 10. — The same as Figu re [SI but for the case of steep mass 
spectrum with a = 4.3 (cf. Fig. I8d). 



effects. Note that the eccentricity dispersion s is virtu- 
ally independent of time above x s h ear for both a = 3.2 
and 4.3 (while inclination dispersion s z slowly grows as 
r inc reas es), which is in complete agreement with equa- 
tion (|C6fl . Finally, numerically computed sq(t) follows 
reasonably well analytical result represented by formula 
lf53*jl . which is demonstrated in Figure 

5. VELOCITY SCALING IN THE PRESENCE OF DISSIPATIVE 
EFFECTS. 

Until now we have been assuming that gravity is the 
only force acting on planetesimals. In reality there should 
be other factors such as gas drag and inelastic collisions 
between planetesimals which might influence their eccen- 
tricities and inclinations. Here we briefly comment on the 
possible changes these effects can give rise to. 

In the presence of gas drag eccentricity and inclination 
of a particular planetesimal with mass to and physical 



size r p tend to decrease with time. Adachi et al. (1976) 
investigated the damping of planetesimal eccentricities 
and inclinations assuming gas drag force proportional to 
the square of planetesimal velocity relative to the gas 
flow. They came out with orbit-averaged prescription for 
the evolution of random velocity of planetesimals which 
can be written at the level of accuracy we are pursuing 
in this study (dropping all possible constant factors and 
assuming that eccentricity is of the same order as incli- 
nation) as 



5e2 
~dt 



Pp r p 



(56) 



where p g is the mass density of nebular gas, p p is the 
physical density of planetesimal, and r\ = (c s /£la) 2 <C 1, 
with c s being the sound speed of nebular gas. In our 
notation (see @ this equation translates into 



ds 2 
dr 



-(s 2 x 



-1/3 



S + J] 



Mc 
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(p P a 3 ) 2 
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1/3 



1 km s 



m 
lQlSg 



1/3 



(57) 



where E g is the surface mass density of the gas disk, 
and numerical estimate is made for p p = 3 g cm~ 3 and 
M c = Mq. For protoplanetary nebula consisting of so- 
lar metallicity gas one would expect E g /S p « 50 — 100, 
but during the late stages of nebula evolution this ratio 
should be greatly reduced by gas removal 7 . 

We use the results of 21 to discuss the effect of the gas 
drag on the planetesimal random velocities for mass dis- 
tributions with a < 3. We know from §4.1 and §4.2 that 
stirring in this case is dominated by biggest planetesimals 
and stirring rate is ~ M<z (r) / s 2 [see equations H28(l , (13811 , 
and l|41|l ]. Comparing this heating rate with damping 
due to gas drag (|5?|) we immediately see that gas drag 
can become important for small planetesimals if planetes- 
imal velocities are large (s 3> 1)- In this case dynamical 
excitation by biggest planetesimals is balanced by gad 
drag. Very interestingly, it turns out that such a balance 
leads to s cx x 1 ^ 12 or x 1 ^ 15 (depending on whether s is 
bigger or smaller than ^(AIc/too) 1 / 3 ), i.e. smaller plan- 
etesimals have smaller random velocities. Such a behav- 
ior has been previously seen in coagulation simulations 
including effects of gas drag on small planetesimals (e.g. 
Wetherill & Stewart 1993). The time dependence of ve- 
locities of small planetesimals affected by the gas drag 
should also be different from that given by l|29|l or (|40p. 
In Figure ITU we schematically display the scaling of s as 
a function of x in the presence of gas drag for 2 < a < 3 
and different values of £. 

Inelastic collisions between planetesimals become im- 
portant for their dynamical evolution roughly when rel- 
ative velocity at infinity with which two planetesimals 
approach each other becomes comparable to the escape 
speed from the biggest of them; in this case gravitational 

7 Note that for planet formation in some exotic environments 
such as early Universe prior to metal enrichment or metal-poor 
globular clusters this ratio can be much larger than 100. 
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Fig. 11. — Schematic representation of the effect of gas drag on 
planetesimal velocity spectrum in the case 2 < a < 3 (cf. Figure 
0. Velocity dependence on planetesimal mass is d ispl ayed for two 
values of gas drag parameter f defined by equation 1571 : £i (dashed 
curve) and £2 (solid curve), with £1 < £2- Velocity spectra overlap 
at high masses for which gas drag is no longer important. Inelastic 
collisions between planetesimals would affect their velocity spec- 
trum in a similar manner. 



focusing is inefficient. The escape speed v esc w 1 m 
s _1 r p /(l km) (for physical density p = 3 g cm -3 ) is much 
larger than the Hill velocity for the same mass flRn ~ 5 
cm s~ 1 a~^J 2 r p /(l km) (for M c = M Q ), which justifies 
initially our previous neglect of inelastic collisions (es- 
pecially far from the Sun where HRh is small). Later 
on, however, continuing dynamical heating in the disk 
would certainly bring velocities of small planetesimals 
above their escape speeds. When this happens, collisions 
with bodies of similar size lead to strong velocity dissi- 
pation. This damping should produce similar effect on 
planetesimal random velocities as gas drag does. Indeed, 
kinetic energy losses in physical collisions between high 
velocity planetesimals are roughly proportional to their 
initial kinetic energies, which is similar to the behavior 
of the gas drag losses. 

At the same time, velocities of light bodies should still 
be smaller than the escape velocity from the biggest plan- 
etesimals, since otherwise these planetesimals would not 
be able to dynamically heat the disk and it would "cool" 
below their escape velocity. Thus, planetesimal veloci- 
ties cannot exceed the escape speed from the surfaces of 
bodies doing most of the stirring, result dating back to 
a classical study by Safronov (1972). This implies that 
gravitational stirring must still be done by the highest 
mass planetesimals (if a < 3) in accordance with what we 
assumed in the case of gas drag. As a result, one would 
again expect s to exhibits a power law dependence on 
mass with positive (but small) slope. Exact value of the 
power law index of this dependence is determined by the 
scaling of energy losses with the mass of planetesimals 



involved in collision 8 . Effect of the collisions should be 
most pronounced for small planetesimals for which rela- 
tive encounter velocities can be much higher than their 
v esc . Threshold mass at which highly inelastic collisions 
become important must constantly increase in time as 
planetesimal disk is heated up by massive planetesimals. 
One should also remember that planetesimal fragmenta- 
tion in high velocity encounters may become important 
in this regime which would significantly complicate sim- 
ple picture we described here. 

In summary, dissipative effects are unlikely to affect 
the dynamics of massive planetesimals. However, they 
become important for small mass planetesimals and lead 
to positive correlation between random velocities of plan- 
etesimals and their masses (unlike the case of pure grav- 
ity studied in which has been previously observed 
in self-consistent coagulation simulations (Wetherill & 
Stewart 1993). It is thus fairly easy to modify our simple 
picture developed in ^to incorporate the effects of gas 
drag and inelastic collisions. 

6. DISCUSSION. 

Although our study has concentrated on a particular 
case of power-law size distributions (most of the mass 
is locked up in the power-law part of the planetesimal 
mass spectrum) its results have a wider range of appli- 
cability. The reason for this is that the velocity scalings 
derived in ^1 depend only on what part of planetesimal 
spectrum dominates stirring and friction and not on the 
exact shape of mass distribution. For example, suppose 
that f(x,r) does not have a self-similar power law form 
of the kind studied in §4.2 (intermediate mass spectrum) 
but its first moment Mi is still dominated by the small- 
est masses (most of the mass is locked up in smallest 
planetesimals, which dominate friction), while the sec- 
ond moment Mi is mainly determined by highest masses 
(heating is determined by biggest planetesimals ~ x c ). 
Then it is easy to see from our discussion in §4.2.1 that 
the velocity distribution has a form of a broken power- 
law in mass, changing its slope from at small masses 
to —1/4 at high masses [see ()4(JII and Q42[l] . If, on the 
contrary, most of the mass is not in smallest bodies but 
in largest ones (Mi and M2 converge near the cutoff of 
planetesimal spectrum, ~ x c ), then both friction and 
stirring are dominated by biggest planetesimals and ve- 
locity spectrum must have zero slope for the entire range 
of planetesimal masses, in accord with our consideration 
of shallow mass distributions (§4.1.1). Whenever both 
stirring and cooling of planetesimals of mass m are dom- 
inated by planetesimals of similar mass, we go back to 
the Case 1 of §4.3. Finally, in all cases when stirring 
and friction are produced by planetesimals much smaller 
than those under consideration, velocity profile scales like 
m^ 1 / 2 (energy equipartition, see Case 2 in §4.3). Simi- 
lar rules can be derived on the basis of the results of jjj] 
for more complicated mass spectra, and this makes our 
approach very versatile. 

One of the interesting features found in our study of 
purely gravitational scattering (see is the develop- 

8 Note that in the crude picture of collision energy losses men- 
tioned above da^/dt is independent of planetesimal mass. Then 
the power law index of velocity spectrum of small planetesimals is 
exactly zero. More accurate treatment of planetesimal collisions 
might result in some mass dependence of random velocities. 
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merit of a plateau at the low mass end of planetesimals 
velocity distribution for shallow and intermediate mass 
spectra (a < 3). These mass spectra are among the 
most often realized in the coagulation simulations which 
makes this prediction quite important for interpretation 
of numerical results. One should be cautious enough not 
to confuse these plateaus with the manifestations of dis- 
sipative processes such as gas drag or inelastic collisions 
which tend to endow velocity distribution with only very 
weak (positive) dependence on mass, see SjSJ 

Another observation which can be made based on our 
results is that the bulk of planetesimals (power law part 
where most of the mass is locked up) rarely reaches a 
complete cquipartition in energy. Equipartition is real- 
ized when a e ^ oc mT 1 / 2 , or in our notation s, s z oc x -1 / 2 . 
Our calculations (§4.3) show that when x < x c such a 
velocity spectrum is only possible for very steep mass 
distributions, a > 4, when both dispersion-dominated 
stirring and friction are due to the smallest planetesi- 
mals with masses m ~ mo (x ~ 1). Mass distributions 
that steep are not often encountered in self-consistent 
coagulation or N-body simulations (Wetherill & Stewart 
1993; Kenyon 2002; Kokubo & Ida 1996, 2000) and yet 
the equipartition argument is very often used to draw im- 
portant conclusions about the bulk of planetesimal pop- 
ulation (see below). 

At the same time, random energy equipartition is ubiq- 
uitous for runaway bodies in the high mass tail (x > x c ) 
if they still interact with the rest of planetesimals in the 
dispersion-dominated regime. Since, by our assumption, 
runaway tail contains very little mass, this equipartition 
is very similar to the case of a > 4, because in both situ- 
ations stirring as well as friction are driven by planetes- 
imals much lighter than the bodies under consideration. 
It is also worth mentioning that this conclusion changes 
when runaway bodies start interacting with small plan- 
etesimals in the shear-dominated regime (see e.g. §4.1.2). 

Although we did not follow the evolution of the plan- 
etesimal mass spectrum self-consistently we can already 
constrain some of the theories explaining mass distribu- 
tions in coagulation cascades. For example, N-body sim- 
ulations of gravitational agglomeration of 10 3 — 10 4 mas- 
sive planetesimals by Kokubo & Ida (1996, 2000) pro- 
duce roughly power law mass spectrum with an index of 
w —2.5 within some range of masses. To explain this 
result Makino et al. (1998) explored planetesimal growth 
in the dispersion-dominated regime with strong gravita- 
tional focusing. They postulated epicyclic energy to be in 
equipartition and found that planetesimal mass distribu- 
tion has a power-law form with a slope of —8/3 w —2.7. 
However, using the results of §4.2.1 one can immediately 
pinpoint the contradiction with the basic assumptions 
which lead to this conclusion. Indeed, slope of —8/3 cor- 
responds to the intermediate mass spectrum (see §4.2) 
which can never reach energy equipartition. Moreover, 
this mass spectrum cannot even have a power-law ve- 
locity distribution with a continuous slope necessary for 
a theory of Makino et al. (1998) to work. As a result, 
the explanation provided by these authors for the mass 
spectrum found in N-body simulations cannot be self- 
consistent. More work in this direction has to be done. 

Although theory developed in 21 is in good agreement 
with our numerical calculations, one should bear in mind 



that our analytical results are asymptotic by construc- 
tion, i.e. they are most accurate when planetesimal disk 
has evolved for a long time and when planetesimal size 
distribution spans a wide range in mass (this can be eas- 
ily seen by comparing Figures El and Q| ■ Real protoplan- 
etary disks during their evolution might not enjoy the 
luxury of having such conditions been fully realized; this 
can complicate quantitative applications of our results 
to real systems. We believe however that this does not 
devalue our analytical results because they provide basic 
understanding of processes involved in planet formation, 
allow one to clarify important trends seen in simulations, 
and enable quick and efficient classification of possible 
evolutionary outcomes for a wide range of parameters of 
protoplanetary disks. 

Several previous studies have concentrated on explor- 
ing velocity equilibria in systems with nonevolving mass 
spectra in which inelastic collisions between planetesi- 
mals acted as damping mechanism necessary to balance 
gravitational stirring (Kaula 1979; Hornung, Pellat, & 
Barge 1985; Stewart k Wetherill 1988). We have out- 
lined a qualitative picture of the effects of dissipative 
processes in jj^l and it agrees rather well with these pre- 
vious investigations. N-body calculations and "particle 
in a box" coagulation simulations represent another class 
of studies in which velocity distributions have been rou- 
tinely computed. As we already commented in §1] these 
calculations usually absorb a lot of diverse physical phe- 
nomena which makes them not easy targets for compar- 
isons. Still, there are several generic features almost all 
simulations agree upon such as (1) roughly power- law 
mass spectra typically with slopes a < 3 within some 
mass range, (2) velocity dispersion is constant or slowly 
increases with mass for small planetesimals, and (3) it de- 
creases with mass for large bodies. This general picture 
is in good agreement with our analytical predictions. 

7. CONCLUSIONS. 

We have carried out an exhaustive census of the dy- 
namical properties of planetesimal disks characterized by 
a variety of mass distributions. Here we briefly summa- 
rize our results for the case velocity evolution driven by 
gravitational scattering of planetesimals. 

Whenever planetesimal mass distribution in a disk has 
a power law slope shallower than —2, we find that plan- 
etesimal velocity is constant up to the upper mass cut- 
off. For slopes between —2 and —3 velocity dispersions 
are constant at small masses but switch to m -1 / 4 de- 
pendence above some intermediate mass. As a result, 
mass spectrum exhibits a pronounced "knee" at this 
mass. Distributions with slopes between —3 and —4 
have purely power-law velocity spectra with slope —f3 
satisfying 1/4 < j3 < 1/2, which can be determined nu- 
merically using equation (|C2(I . Finally, planetesimals in 
disks characterized by mass spectra steeper than m~ 4 
are in energy equipartition, i.e. velocity dispersions scale 
as m -1 / 2 . Such a difference in behaviors is caused by 
the fact that gravitational stirring and dynamical fric- 
tion receive major contributions from different parts of 
planetesimal mass spectra in different cases. 

We also consider a possibility that mass distribution 
has a tail of massive "runaway" bodies sticking out be- 
yond the exponential cutoff of the power-law mass spec- 
trum. We have found that the low mass end of the tail 
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experiencing dispersion-dominated scattering by all plan- 
etesimals exhibits complete energy equipartition: a e .i oc 
mT x l 2 . Most massive runaway bodies interact with all 
planetesimals in the shear-dominated regime, which leads 
to highly anisotropic velocity distributions of these bod- 



ies: a e <x m x / 6 while cr. 



oc m 



-1/2 



Time dependence of 



velocities is determined mainly by the evolution of the 
cutoff of the power law mass spectrum moX c (r). These 
asymptotic results derived using analytical means are in 
good agreement with numerical calculations of planetes- 
imal dynamical evolution presented in 2J They can be 
easily generalized to cover more complicated planetesi- 
mal mass distributions. 

Finally, we investigate qualitatively the impact of 
damping processes such as gas drag or inelastic collisions 
on the planetesimal disk dynamics. We find that manifes- 
tations of these effects are only important for small mass 
planetesimals; they exhibit themselves in the form of gen- 



tle increase of planetesimal random velocities with mass 
which has been previously observed in "particle in a box" 
coagulation simulations (Wetherill & Stewart 1993). 

Future work in this direction should target the self- 
consistent coupling of the evolution of planetesimal ve- 
locities to the evolution of mass spectrum of planetesi- 
mals due to their coagulation. Successful solution of this 
problem would greatly contribute to our understanding 
of how terrestrial planets grow out of swarms of plan- 
etesimals in protoplanetary disks. 
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APPENDIX 



NUMERICAL PROCEDURE. 



Evolution of planetesimal eccentricities and inclinations is performed by numerically evolving in time the full system 
with scattering coefficients smoothly interpolated between shear- and dispersion-dominated regimes. Planetesimals 
of different masses are distributed in logarithmically spaced mass bins, with bodies in each bin being more massive than 
planetesimals in previous bin by 1.2. Using numerical orbit integrations we have determined the values of constant 
coefficients in @ relevant for the shear-dominated scattering to be 



B x fa 11.4, B 2 ps 4.0, D x « 11.4, D 2 ~ 4.2. 



(Al) 



Scattering coefficients in the dispersion-dominated regime are taken from Stewart & Ida (2000) who derived analytical 
expressions for viscous stirring and dynamical friction rates in the two-body approximation. We translate their results 
for our gravitational stirring and friction functions H\ 2 and K\ 2 . This provides us with the dependence of coefficients 
Ai <2 and C\ i2 in © on the inclination to eccentricity ratio <Ji/o e . Constants a% t 2 in (J7|) are fixed to be a\ = 1.7, a 2 = 1. 

We assume that mass scale evolves as x c (t) = F + (er)* where F ps 40, e = 10~ 4 , and X — 4 in the case of 
shallow mass spectrum with a = 1.5 and x — 1-5 in all other cases. Small parameter e is introduced to mimic the 
difference between the timescale of dynamical evolution of planetesimal disk and the timescale of planetesimal growth 
(which is usually much longer). Besides this, such a form of x c (r) is chosen arbitrarily. Runaway tail is assumed to 
have a power-law form with the same slope a but much smaller normalization than planetesimal mass distribution. 
Tail always extends 8 orders of magnitudes beyond x c to allow interesting velocity regimes to fully develop. Initial 
distribution of planetesimal eccentricities and inclinations is set to s(cc,0) = s z (x,0) = 3x~ 1//2 . 

VELOCITY EVOLUTION IN A "MIXED" STATE FOR INTERMEDIATE MASS SPECTRUM. 

In the case of intermediate mass spectrum runaway bodies with masses x between Xh [defined in (|44|l ] and x s h e ar 
interact with small planetesimals [less massive than x s (x) defined by (|45|) ] in the dispersion-dominated regime, but 
with large ones (heavier than x s (x)) in the shear-dominated regime. As a result, the equations for s and s z are now 
hybrid versions of l|30fl and (|35[1 : 



s (x) 



, + x* 2 f(x*) „ „ 2 f , + x*f(x*) C 
dx* J -±-^- - 2A 2 s 2 x / dx* 17 ; 4 ' 



i 

x s (x) 



x s (x) 

I" 

1 

x s (x) 



x 



2/ 



/ dx* X -^£l-2B 2 slx J dx* X -^fl + ^L J ^VV(^) , 

1 1 x s (x) 



^ I dx*x* 2 f(x*)~^- J dx*x*f(x*), (Bl) 

dx*x* /(x*XB2) 



x s (x) 



x s (x) 

2D 2 s 2 



s (x) 



Following the approach taken in §4.1.2 we neglect l.h.s. of both equations and use (|4U|I and 142f) to evaluate integrals 
in evolution equations. After careful comparison of different contributions in the r.h.s. we find that 



C 1 ^ps2A 2 M 1 



2 

s z x 



B 



2-2/3 

7/2-a 

2 1/2 
s Q x k 



+ Di 



2 !/ 2 
3-4/3 



OC 

/ dx*(x*) s / 2 f(x*)n2B 2 M 1 ^. 

J 4 



(B3) 
(B4) 



18 



Equation l|B3(l implies that eccentricity stirring of runaway bodies in the mass range ^ < i< x s hear is done mainly 
by planetesimals with masses ~ x c , for which the interaction proceeds in the s/iear-dominated regime. This heating is 
balanced by dynamical friction due to the smallest planetesimals, lighter than Xk, which interact with these massive 
bodies in the dispersion-dominated regime. 

Vertical heating is somewhat different. Gravitational friction is again dominated by smallest planetesimals which 
are in the dispersion-domin ated mode. Stirring, however, depends on shape of the mass distribution. If a < 5/2, 
second term in the l.h.s. of (|B4(1 — shear-dominated heating by planetesimals with masses between ~ x s and ~ x c 
— dominates 9 over the first term which represents the effect of dispersion- dom inated heating by planetesimals with 
masses ~ x s . I n th e opposite case, when a > 5/2 both terms in the l.h.s. of l|B4|l produce roughly equal contributions. 
Using (|B3ll & i|B4|l one can easily derive expressions (|46f) - (l48ll . 

DETAILS OF VELOCITY EVOLUTION FOR THE CASE OF STEEP MASS SPECTRUM. 

Substituting our guess (5U into © and going through all possibilities appropriate for a > 3 we find that self- 
consistent solutions of power-law type can exist in only two cases: 

• Case 1: 3 - a + 2/3 > & 2-a + 4/?>0. 

• Case 2: 3 - a + 2/3 < & 2 - a + 4/3 < 0, 
We now consider separately these two possibilities. 

Planetesimal velocities. 

Case 1. 

In this case one can easily see that integrals in the r.h.s. of (|26|l are mostly contributed by x* ~ x. This allows us 
to rewrite l)26fl in the following form: 



x _ 20 dsl _ x*-«+W [df-^ + V 



OO 

t , 1 



dr s 2 J l + t-W 



o 



Ai - 2——r 7r A 2 



1 + t 1 1 + t-W 



(CI) 



In arriving at this expression we have used the fact that f(x, r) = x~ a for a > 3 and x c , and extended the 

integration range over t = x* jx f rom to oo (because only local region t ~ 1 matters). From this equation we can 
readily see that the l.h.s . of (|C1(I can be neglected for x 3> 1. As a result, a self-consistent power-law solution for s 
exists only if integral in (|C1(I is equal to zero. This leads to the following constraint on the required value of (3: 

h{a,P)= / dt——— M -K, I 2 (a,(3) = / dt „ . (C2) 



I 2 (a,(3) 2A 2 ' ^^'^"7^(1 + ^2/3)2' ^v-.^-y - 1 + t -2 - 



Solving this equation for given a a nd A i/A 2 one can find the slope of the mass spectrum (3(a, Ai/A 2 ). At a first sight 
it is not at all clear that equation (|C2|) should in general possess a solution for (3. However, closer look at the problem 
reveals some interesting patterns. 

First of all, it follows from (JSJ) that r.h.s. of (|C2|I has to be bigger than 1. Second, in the case of energy equipartition 
Ii(a , l/2)/I 2 (a, 1/2) = 1. Third, Ii(a,{3) — > oo as (3 — > (a — 2)/4, while I 2 stays finite. These observations prove that 
(|C2|I always has a solution for (3 satisfying the constraint posed by equation (|55|l . Combining restriction (|55|l with 
initial constraints of Case 1 we find that Case 1 is only possible for mass spectra with power law slope 3 < a < 4. 

In fact, one can do things even better by considering separately eccentricity and inclination scalings with mass, i.e. 
using both equations i|16|) . Assu ming that the ratio of inclination to eccentricity s z /s is still constant for x < x c , on e 
would obtain in addition to l|C2|) another equation dictated by the inclination evolution. It would be identical to (|C2(1 . 
but with Ai t2 replaced by B\_ 2 . Since both A\ >2 and B\_ 2 depend only on s z /s (see @ one would have two equations 
for two unknowns: (3 and s z /s. Solving them one can uniquely fix both the power law index of planetesimal velocity 
dependence on mass and the ratio of inclination to eccentricity of planetesimals (which is left undetermined in our 
simplified analysis of 

Case 2. 

Restrictions imposed on a and (3 by the conditions of Case 2 imply that all integrals in the r.h.s. of (|26|l are 
dominated by the lower end of their integration range, i.e. by masses ~ 1. Neglecting time derivative in (|26Jl and 
balancing contributions in the r.h.s. we find that (3—1/2 and 

s(x, t) ~ sq(t)x~ 1 / 2 , X < x c - (C3) 

This solution demonstrates that for x 3> 1 our omission of the l.h.s. of (|26|l is justified for arbitrary behavior of x c (t). 
It also imposes important restriction on the mass spectra for which this scaling can be realized: a mus be bigger than 
4. This means that we have found solutions for all possible positive values of a and our study is complete in this sense. 

9 For a < 5/2 integral in IB4I converges at the upper end of its range, near x c . When a = 5/2 this in tegra l is contributed roughly 
equally by equal logarithmic intervals in mass between ~ x s and ~ x c ; then the second term in the l.h.s. of IB4I dominates over the first 
one by ~ ln(x c /x s ). 
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Velocities of runaway bodies. 

Case 1 (3 < a < 4). 

Velocity evolution of runaway bodies in the Case 1 is similar to the situation with the intermediate mass spectrum. 
One finds that when mass x satisfies x c < x < x^ = SqX~ 3 @ velocity profile is shaped by the dispersion-dominated 
interaction with all small mass planetesimals (x < x c ): 

„( x ,T)K8,(x,T)a^(xl-w¥^] , 3<a<4, x c < x < x h . (C4) 

The biggest contribution to both stirring and dynamical friction comes the most massive planetesimals with masses 
^ x c . 

Going to heavier bodies, x^ < x < x s h ea r = Sq, one finds scattering to be in mixed state: some planetesimals [those 
heavier than x s (x) = (sq/x 1 ^ 3 ) 1 ^} interact with big bodies in the shear-dominated regime while the other part (those 
lighter than x s ) is in the dispersion-dominated regime relative to runaway bodies. It turns out that when 3 < a < 4 
the biggest contributors to both the heating and friction of big bodies are planetesimals of mass ~ x s . Since this is 
just at the boundary bet ween the shear- and dispersion-dominated regimes, one can expect s z to behave in the same 
way as s does [similar to (|C4(l ], Indeed, we find this to be the case (Figure |HK) : 

/ / n s V(2/3) 

s(x, t) ps s z (x, t) ps s /6 ps ^ +1//3 y 6 , 3 < a < 4, x h < x < x shear . (C5) 

Runaway bodies with x > x s h ear = Sq ps r 3 / 4 experience shear-dominated scattering by smallest planetesimals; 
using (|36Jl one finds that 

s(x,t) ps x~ 1/6 I —j , s z (x,t) ^ s (t)x~ 1/2 , a > 4, a; > x shear . (C6) 

Case 2 (a > 4). 

We saw in §C.l that in the Case 2 planetesimal velocity evolution is determined purely by the smallest bodies. As a 
result, it does not matter whether one studies the dynamics of planetesimals or runaway bodies — all the conclusions 
of §C.l stay unchanged as long as the interaction with smallest planetesimals occurs in the dispersion-dominated 
regime. Thus we predict that for x < x s fi ear velocity spectrum is still given by (|C3|) . Runaway bodies of l arger 
mass, x > x s } lear , feel shear-dominated scattering by all planetesimals and their velocity dispersions are given by lC6t 
(exactly like in Case 1). 
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Table CI. Summary of analytical results for the velocity scaling 

WITH MASS. 



Range of a Mass interval 
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